Spin-orbit coupled weakly interacting Bose-Einstein condensates in harmonic traps 
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We investigate theoretically the phase diagram ol a spin-orbit coupled Bose gas in two-dimensional 
harmonic traps. We show that at strong spin-orbit coupling the single-particle spectrum decomposes 
into different manifolds separated by h\u±, where lu± is the trapping frequency. For a weakly 
interacting gas, quantum states with skyrmion lattice patterns emerge spontaneously and preserve 
either parity symmetry or combined parity-time-reversal symmetry. These phases can be readily 
observed in a spin-orbit coupled gas of 87 Rb atoms in a highly oblate trap. 
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Spin-orbit (SO) coupling leads to many fundamental 
phenomena in a wide range of quantum systems from 
nuclear physics, condensed matter physics to atomic 
physics. For instance, in electronic condensed matter sys- 
tems SO coupling can lead to quantum spin Hall states 
or topological insulators [J, which have potential ap- 
plications in quantum devices. Recently, SO coupling 
has been induced in ultracold spinor Bose gases of 87 Rb 
atoms by the so-called "synthetic non-Abelian gauge 
fields". Combined with unprecedented controllability of 
interactions and geometry in ultracold atoms, this manip- 
ulation of SO coupling opens an entirely new paradigm 
for studying strong correlations of quantum many-body 
systems under non-Abelian gauge fields. 

In this context, over the past few years there have been 
great theoretical efforts to determine quantum states of 
an SO coupled spinor Bose-Einstein condensate (BEC) 
In a recent work by Wang et al. two distinct 
phases are identified for a homogeneous two-dimensional 
(2D) spin-1/2 BEC. Depending on the relative magnitude 
of intra-species (g) and inter-species (<? n ) interactions, all 
bosons can condense into either a single plane- wave state 
(<7 < g n ) or a density-stripe state [g > <7 n ). 

The purpose of this Letter is to show that the presence 
of a harmonic trap, which is necessary in experiments, 
can change dramatically the phase diagram of SO coupled 
BECs. At strong SO coupling the single-particle spec- 
trum decomposes into discrete manifolds, analogous to 
discrete Landau levels. Non-trivial quantum states with 
skyrmion lattices emerge when all bosons occupy into the 
lowest manifold (LM). These properties are fundamen- 
tally different from that of a homogeneous system. We 
note that, in a previous work, the NIST group has exper- 
imentally realized an artificial Abelian gauge field which 
leads to the observation of vortex lattice in a non-rotating 
87 Rb condensate [loj . Our work represents an important 
extention into the regime of non-Abelian gauge field in 
which the spin degrees of freedom play an essential role. 

Our main results are summarized in Fig. 1, which 
shows the ground state as functions of interatomic inter- 
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Figure 1: (color online), (a) Phase diagram of a trapped 2D 
BEC with a strong SO coupling A = 20, where the single- 
particle spectrum forms discrete manifolds. For the weak 
interaction considered here, only the LM is occupied. The 
phases I and II preserve, respectively, the parity and parity- 
time-reversal symmetries. There are several sub-phases in- 
dicated by A, A' and B, which differ in the density profile 
and/or angular momentum. The mean-field density patterns 
in different phases of spin-up bosons are shown in Fig. El (b) 
and (c) Phase diagram at weak SO coupling. Here the phases 
are determined without the restriction to the LM approxima- 
tion. The insets illustrate the density profiles of the two spin 
components in phases IA and IIA. 



action at a dimensionless SO coupling strength A. By us- 
ing mean-field theory and exact diagonalization, we find 
that: (i) The ground state falls into two classes of quan- 
tum phases, I and II, preserving respectively the par- 
ity (V) and parity-time-reversal (VT) symmetries. Both 
symmetries are satisfied by the model Hamiltonian [see 
Eqs. (JT|) below], (ii) In each class, there are several sub- 
phases (IA, IA', IB and IIA, IIB) differing in the den- 
sity distribution and/or total angular momentum, (iii) 
The transition between different phases depends on inter- 
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atomic interactions. At weak intra-species interactions 
below a critical value, g < g c , the ground state is a half- 
quantum vortex state (IA) if g < g n and a superposi- 
tion of two degenerate half-quantum vortex states (IIA) 
otherwise. The phases I A and IIA vanish in the limit 
of strong SO coupling, but dominate the phase diagram 
in the opposite. When the intra-species interactions be- 
comes larger (<? > g c ), there is an interesting reverse of 
the symmetry class, i.e., interactions change the phase 
I A into IIB and the phase IIA into I A' and then IB. In 
the phases IIB and IB, skyrmion lattices emerge sponta- 
neously without rotation, (iv) At g = g.,, the phases are 
ordered by quantum fluctuations. Using exact diagonal- 
ization, we find that the phases follow those at g < g u . 

Model Hamiltonian and energy spectrum. - We con- 
sider iV-bosons in a 2D harmonic trap V(p) = MuP\_p 2 /2 
with a Rashba SO coupling V so = —iXu(d y a x — d x a y ), 
where cj x ^ y ^ z are the Pauli matrices. The model Hamilto- 
nian is given by % = Jio + Hint , where 

H a = J dr^+ [-h 2 V 2 /(2M) + V(p) + V so ] *,(la) 

Hint = fdr \{g + g n )h 2 + (g - g u )S 2 } /4 , (lb) 



states of VT- 



\& = [^(r), ^(r)] 7 " denotes collectively the spinor Bose 
field operators, and n,S z = *+vfi t ±\]/+\]/ i . We define 

two characteristic lengths, a± = y/H/(Muj±) for the har- 
monic trap and a\ — ft 2 / (MXr) for the SO coupling. The 
dimensionless SO coupling strength can be then defined 
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The Hamiltonian 



as A = a^/ax = (M/H 3 )^ 2 X r /lu 
is invariant under two symmetry operations, associated 
respectively with the anti-unitary time-reversal operator 
T = icr y C, where C takes the complex conjugate, and the 
unitary parity operator V = a z T, where X is the spatial 
inversion operator. The Hamiltonian is also invariant un- 
der the combined VT operator, which is unitary since V 
and T anti-commute with each other, i.e., {V,T} = 0. 

In polar coordinates (p, ip) , the single-particle eigen- 
wavefunctions of may be written in the form, $ m (r) = 
[(/> t (p)e m ' p , (/)i(p)e^ m+1 ^] T , which is energetically de- 
generate with its time reversed partner T$ m (r) = 
l(t>l(p)e- %( - m+1 ^,-<j) t (p)e- tmv } T . This degeneracy is a 
direct consequence of the Kramers' Theorem. Here we 
may restrict m to be non-negative integers, as a nega- 
tive m state can be regarded as the time reversal part- 
ner for a state with m > 0. In this construction, <E> m 
and 7~&m are both parity eigenstates with corresponding 
eigenvalues (—1)'™ and (— l) m+1 , respectively. However, 
they break the VT symmetry. The lowest single-particle 
state occurs at to = and has a half-quantum vortex 
configuration j^. Due to the degeneracy, any linear su- 
perposition of <I> m and T&m — which breaks the parity 
symmetry — is also an eigenstate of the system. In par- 
ticular, we may choose the equal-weight superposition as 
(&m + T$ m )/\/2 which can be easily shown to be eigen- 




Figure 2: (color online), (a) Single-particle energy spectrum. 
The lines show the empirical Eq. ((2)1 . (b) The W-mnction for 
the lowest four single-particle states in the LM. 

The wavefunctions and the corresponding eigenener- 
gies can be found numerically. At large SO coupling 
(i.e., A > 5), to a good approximation we find numer- 
ically that the low-lying spectrum forms discrete mani- 
folds with spacing hiu± (indexed by an integer n > 0), 



-A 2 + (2n + l)+m(m + l)/A 2 Hlu ± /2 . (2) 



There are about 2y/2 X levels within each manifold with 
the smallest level spacing AE = huj±/X 2 . The discrete 
manifolds of spectrum are similar to the well-known Lan- 
dau levels, formed when a charged particle moves in mag- 
netic fields. However, the reasons for their formation are 
very different. In our case of large SO coupling, without 
trap the spectrum is characterized by a continuous mo- 
mentum k and is given by ek = [— A 2 /2+ (fc± X) 2 /2]haj±, 
with infinite degeneracy along the azimuthal direction. 
The inclusion of trapping potential quantizes the ra- 
dial motion for k and the azimuthal motion, giving the 
standard quantization contribution of (n + l/2)fou)± and 
(to + 1/2)" /(2X 2 )haj± to the energy, respectively. 

For a weakly interacting BEC with gN, g^N <C fiio±, 
only the LM is occupied. It is thus convenient to expand 
the field operator '5 = J2 m ^m(r)a m , where < I ) m (r) is the 
single-particle wavefunctions at the LM with energy e m . 
The many-body Hamiltonian may then be rewritten as, 



rn (l rn Cl rn 



^ Vijkiaj aj a A; a/, 

ijkl 



(3) 



where the interaction elements Vy/y can be calculated 
straightforwardly for the contact interatomic interac- 
tions. We solve Eq. ([3]) numerically by using both mean- 
field theory llj and exact diagonalization [12j , for a con- 
served total angular momentum X)m( m + ^/^) a m a m = 
Nmtot- Within mean-field, we replace a m by a complex 
number N x l 2 c m and minimize the GP energy Eqp/N = 
E m e m | Cm I + {N-1) J^ijki Vijkic*c*c k ci, under the con- 



straints J2 m \ c r, 



1 and E m (™+V2) \c n 



m tot . 
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In practice, we truncate the angular momentum to |m| < 
m c (up to m c = 16). 

Symmetry of condensate states. - In the presence of 
the interaction represented by Eq. (|lb[) , the many-body 
Hamiltonian still possesses both V and VT symmetries. 
As we have shown above, for a non-interacting system, 
we may choose the single-particle ground state to be an 
eigenstate of V, or of VT, or of neither operator. In the 
mean-field level, this freedom of choosing different sym- 
metry eigenstates may be removed by inter-atomic in- 
teractions. In other words, the symmetry of condensate 
states would be determined spontaneously by interaction. 
We have found that in the weakly interacting limit we are 
interested in here, the ground state is either an eigen- 
state of V, or that of VT- Which symmetry the ground 
state will possess can be determined in the following way. 
Let us consider an eigenstate of V with wavefunction 
^■p = [cj>^(r), cf>i(r)] T . The corresponding eigenstate of 
VT can be constructed as $707- = (§ v ±T$v)/V2- The 
mean-field energy difference between these two states is 
determined by the S 2 term in Eq. (|lb[) which breaks the 
spin rotational symmetry in the interaction Hamiltonian: 

AE sp ($) = E($ VT ) - E($-p) = (ffu - 9)W (*) /4, (4) 

where W (*) = / dr[( |0 t | 2 - I0i| 2 ) 2 - (Mi + Wl- 
The ground state will be a P-eigenstate if AE sp ($) > 
for which we have n a {v) = n a (—r), or a PT-eigenstate 
if AE sp ($) < for which we have n-f(r) = r). The 
VK-functions of several parity eigenstates are shown in 
Fig. Efb). Equation ([J} also shows that the symmetry of 
the ground state is sensitive to the relative magnitude of 
the interaction parameters g and g^. 

Phase diagram in the LM. - Our symmetry argument 
suggests that all the condensate states could be classified 
by its V or VT symmetry, to be referred to respectively 
as phases I and II hereafter. We now check numerically 
this argument in the quantum Hall like regime with all 
bosons occupying into the LM, as shown in Fig. [lja) for 
A = 20. The characteristic density distributions for spin- 
up bosons in each phase are shown in Fig. [3] 

At sufficiently weak interactions, where the character- 
istic interaction energy g(N — l)a]_ is smaller than the 
lowest intra-manifold spacing AE = hu±/X 2 , only the 
ground single-particle state is occupied. The condensate 
state is thus either half-quantum vortex states of $0 ( or 
T$o) or their superposition. As W($o) > as shown 
in Fig. [^b) , we conclude that the ground state is a VT- 
eigenstate for g > g^ (HA) and it is a half-quantum vor- 
tex state (a 'P-eigenstate) for g < g^ (IA) . Their spin- up 
density patterns are shown in Figs.^a) and (d), respec- 
tively. 

When the interaction becomes larger, more and more 
single-particle states are occupied. The occupation of 
the first excited single-particle state ($1 and T$i) oc- 
curs at g c (N — l)a\ ~ Q.0367fuj±, where the critical 
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Figure 3: (color online). Density patterns of spin-up bosons in 
the different ground states at three intra-species interactions 
g(N- l)a\: (a,d) 0.02/jwx, (b,e) O.lfiwx, and (c,f) 0.2hu±. 



interaction strength g c is determined from the equation 
c + (JV-l)Vbooo = ei + (^-l) Vim. As W($ m ) <0for 
m > 1, we find an interesting reverse of the phase dia- 
gram when g > g c : the P-preserving phase (I A) changes 
into a PT-preserving phase (IIB) at g < g^, while the 
PT-preserving phase (IIA) changes into a P-preserving 
phase (IA' and IB) if g > Q.2g^. The phases IA' and 
IB differ in the total angular momentum mtot and den- 
sity distribution. In Phase IB, mtot is suppressed to 
zero by large interatomic interactions. Note that in the 
phases (IIB) and (IB), we observe regular lattice pat- 
terns. In particular, a hexagonal lattice form gradually 
in the phase IIB, as shown clearly in Figs. (3Je) and (f). 
In Fig. [31 we show the corresponding spin texture of the 
state, from which one can see that the system represents 
a lattice of skyrmions. Skyrmion lattice can be generated 
by rotating a spinor condensate (l3| . Here the skyrmion 
texture is induced by the SO coupling without rotation. 

The symmetry of the ground state at g = g^ can not 
be determined within mean-field theory, since in this case 
AE sp ( < f ) ) = [see Eq. (g])] and the energy becomes in- 
variant for different m^t. However, it can be ordered by 
quantum fluctuations [4| , which are well captured by ex- 
act diagonalization. We have calculated the energy as a 
function of mtot at g(N — l)a\/ (hjj±) = 0.02 and 0.1 for 
N = 4, 8, and 12. With increasing N, the exact diagonal- 
ization result approaches the mean-field prediction. We 
find that the ground state at g(N — l)aj_ = 0.02haj± has 
a spontaneous angular momentum m to t = — 1/2 or +1/2, 
while the ground state at g(N — \)a\ = 0.17kjj_ occurs 
at mtot = 0. Therefore, we identify that the phases at 
g = follow those at g < g^. This is in agreement 
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Figure 4: (color online). Spin texture S = (l/2)$ + cr$ cor- 
responding to the state represented in Fig. [3jf). The ar- 
rows represent the transverse spin vector (S x , S y ) with color 
and length representing the orientation and the magnitude of 
the transverse spin. The contour plot shows the axial spin 
S, = (1/2) (<$-<%). 



with the result of Ref. jij, which employs a different "or- 
der from disorder" argument. 

Phase diagram beyond LM. - So far we have clarified 
the phase diagram at a particular SO coupling A = 20 in 
the weakly-interacting LM regime. However, the qualita- 
tive picture of diagram may persist beyond the regime of 
LM, as far as our symmetry argument holds. To check 
this, we performed a direct numerical calculation based 
on the full Gross-Pitaevskii (GP) equation derived from 
Eqs. (JlJ without making the LM assumption. In the 
regime as shown in Fig. QJa), the results are in good 
agreement with the LM calculation. At larger interac- 
tion strength when higher manifolds get mixed in the 
ground state, we have found from the GP calculation 
that Phase IIB in Fig.QJa) will change to a density-stripe 
phase with V symmetry, while Phase IB will change to 
a plane-wave phase with VT symmetry. The density- 
stripe and the plane-wave phases have been shown to be 
the mean-field ground state for a homogeneous system 
For the trapped system as studied here, at large 
interaction strength, the effect of the trap becomes less 
important and our results are therefore consistent with 
those reported in Ref. With decreasing A, we antic- 
ipate that the phases IA and IIA will gradually become 
dominant in the diagram, as we find numerically that 
g c oc 1/A 2 increases very rapidly. The skyrmion lattice 
phase, related to the LM formation, may disappear. This 
is confirmed by the GP calculation for smaller SO cou- 
pling and the results are represented in Fig. Q] (b) and 
(c) . The half-quantum vortex state and its superposition 
dominate over a much larger parameter space as com- 
pared to the large SO coupling case. A more detailed 
study of the complete phase diagram and the prop erties 
of different phases will be presented elsewhere 



Experimental relevance. - We finally consider the ex- 
perimental feasibility. A Rashba SO coupling can be in- 
duced in spinor 87 Rb gases ■ The interaction strengths 
of 87 Rb atoms may be tuned by properly choosing the 
parameters of the laser fields that induce the SO cou- 
pling The two-dimensionality in such system has now 
been routinely realized by imposing a strong harmonic 
confinement V(z) = Mu>^z 2 /2 along the z-direction with 
lu z 3> The critical temperature for an ideal 2D SO 
BEC is given by T c = (c\ / 'tt)V '3NHuj ± j 'ks , where the 
prefactor c\ < 1 takes into accout the suppression due 
to the SO coupling. Taking parameters from a recent ex- 



periment [15| with lo± = 2ir x 20.6 Hz and N ~ 10 5 , we 



find at A = 10, c\ ~ 0.6 and kgT c ~ 120 nK. Experimen- 
tally, BEC temperature below 0.5 nK has been recorded 
[lij . which is also lower than huj^/ks- The mean-field 
LLL regime is therefore readily attainable with current 
technologies. 

Conclusion. - In summary, we have investigated the 
phase diagram of a spin-orbit coupled spinor BEC in 
harmonic traps, by using mean-field theory and exact 
diagonalization method. We have predicted that the con- 
densate states preserve the parity or parity-time-reversal 
symmetry and exhibit spontaneous vortex and skyrmion 
lattice structure in the lowest energy manifold which is 
induced by strong spin-orbit coupling. Our results are 
valid for weak correlations with large number of bosons. 
Strongly correlated states, analogous to the fractional 
quantum Hall states, would emerge with small number 
of bosons [l3| ■ These can be addressed using exact diag- 
onalization method in future studies. 
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Note added. - When our manuscript was under re- 
view, we became aware of a preprint [18J , in which the 
authors addressed the same problem at g = g^. 
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